#ifndef SPLINE_FUNCTIONS_H
#define SPLINE_FUNCTIONS_H

#include <vector>
#include <functional>
#include "Eigen/SparseCore"
#include "Eigen/SparseLU"
#include <stdexcept>
#include <cassert>
#include <iostream>
#include <algorithm>

// 基础样条函数类(ok)
class SplineFunction {
public:
    virtual ~SplineFunction() {}
    virtual double evaluate(double x) const = 0; // 评估样条函数
};

// pp形式的样条函数基类(ok)
class PPFormSpline : public SplineFunction {
protected:
    std::vector<double> knots; // 节点
    std::vector<double> coefficients; // 系数

public:
    PPFormSpline(const std::vector<double>& knots, const std::vector<double>& coefficients)
        : knots(knots), coefficients(coefficients) {}

    // 可以添加一些pp形式样条特有的成员函数
};

// B样条函数基类(ok)
class BSpline : public SplineFunction {
protected:
    int degree; // 阶数
    std::vector<double> knots; // 节点
    std::vector<double> coefficients; // 系数

public:
    BSpline(int degree, const std::vector<double>& knots, const std::vector<double>& coefficients)
        : degree(degree), knots(knots), coefficients(coefficients) {}
    virtual double evaluate(double x) const override{
        return 0;
    };
};

// pp形式的线性样条函数类(ok)
class LinearSpline_p : public PPFormSpline {
public:
    LinearSpline_p(const std::vector<double>& knots, const std::vector<double>& coefficients)
        : PPFormSpline(knots, coefficients) {}

    // 重写评估函数
    virtual double evaluate(double x) const override{
        assert(knots.size() > 0 && coefficients.size() == knots.size());
        int N = knots.size();
        for (int i = 0; i < N - 1; ++i) {
            if (x >= knots[i] && x < knots[i + 1]) {
                double slope = (coefficients[i + 1] - coefficients[i]) / (knots[i + 1] - knots[i]);
                return coefficients[i] + slope * (x - knots[i]);
            }
        }
        throw std::out_of_range("x is out of the spline's domain.");
    }
};

// pp形式的周期三次样条(ok)
class PeriodicCubicSpline_p : public PPFormSpline {
public:
    PeriodicCubicSpline_p(const std::vector<double>& knots, const std::vector<double>& coefficients)
        : PPFormSpline(knots, coefficients) {}
        
    // 重写评估函数
    virtual double evaluate(double x) const override{
        assert(knots.size() > 1 && coefficients.size() == knots.size());
    
        int N = knots.size();
        Eigen::VectorXd M(N);
        Eigen::VectorXd lambda(N);
        Eigen::VectorXd mu(N);
        Eigen::VectorXd L(N);

        // 计算 mu 和 lambda 
        for (int i = 1; i < N-1; ++i) {
            double h_i = knots[i] - knots[i-1];
            double h_next = knots[i + 1] - knots[i];
            mu[i] = h_i / (h_i + h_next);
            lambda[i] = h_next / (h_i + h_next);
            if(i == N-2){
                L[i] = (coefficients[0] - coefficients[i]) / (h_next*((h_i + h_next))) - (coefficients[i] - coefficients[i-1]) / (h_i*((h_i + h_next)));
            }
            else{
                L[i] = (coefficients[i+1] - coefficients[i]) / (h_next*((h_i + h_next))) - (coefficients[i] - coefficients[i-1]) / (h_i*((h_i + h_next)));
            }
        }

        // 设置三对角矩阵 A 和向量 b
        Eigen::SparseMatrix<double> A(N, N);
        Eigen::VectorXd b(N);

        // 填充三对角矩阵 A 和向量 b
        for (int i = 1; i < N - 1; ++i) {
            A.insert(i, i) = 2;
            A.insert(i, i - 1) = mu[i];
            A.insert(i, i + 1) = lambda[i];
            b[i] = 6*L[i];
        }
        // 第一个方程的系数
        A.insert(0, 0) = 1;
        A.insert(0, N-1) = -1;
        b[0] = 0;

        // 最后一个方程的系数
        A.insert(N-1, 0) = 2*(knots[1]-knots[0]); 
        A.insert(N-1, 1) = knots[1]-knots[0]; 
        A.insert(N-1, N-2) = knots[N-1]-knots[N-2]; 
        A.insert(N-1, N-1) = 2*(knots[N-1]-knots[N-2]); 
        b[N-1] = 6*(coefficients[1]-coefficients[0])/(knots[1]-knots[0]) - 6*(coefficients[N-2]-coefficients[0])/(knots[N-2]-knots[N-1]);

        // 使用 Eigen 解线性方程组 A * m = b
        Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
        solver.compute(A);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to decompose the matrix.");
        }
        M = solver.solve(b);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to solve the system.");
        }

        for (int i = 0; i < N - 1; ++i) {
            if (x >= knots[i] && x < knots[i + 1]) {
                double h_i = knots[i+1] - knots[i]; 
                double K_i;
                if(i == N-2){
                    K_i = (coefficients[0] - coefficients[i])/h_i;
                }
                else{
                    K_i = (coefficients[i+1] - coefficients[i])/h_i;
                }
                double m_1 = K_i - (M[i+1]+2*M[i])*h_i/6;
                double m_2 = K_i + (M[i]+2*M[i+1])*h_i/6;
                double sx = coefficients[i] + (x-knots[i])*m_1 + (x-knots[i])*(x-knots[i])*(K_i-m_1)/(knots[i+1]-knots[i]) + (x-knots[i])*(x-knots[i])*(x-knots[i+1])*(m_1+m_2-2*K_i)/((knots[i+1]-knots[i])*(knots[i+1]-knots[i]));
                return sx;
            }
        }
        throw std::out_of_range("x is out of the spline's domain.");
    }
};

// pp形式的自然三次样条(ok)
class NaturalSpline_p : public PPFormSpline {
public:
    NaturalSpline_p(const std::vector<double>& knots, const std::vector<double>& coefficients)
        : PPFormSpline(knots, coefficients) {}

    // 重写评估函数
    virtual double evaluate(double x) const override{
        assert(knots.size() > 0 && coefficients.size() == knots.size());
    
        int N = knots.size();
        Eigen::VectorXd M(N);
        Eigen::VectorXd M0(N-2);
        Eigen::VectorXd lambda(N);
        Eigen::VectorXd mu(N);
        Eigen::VectorXd L(N);

        // 计算 mu 和 lambda
        M[0] = 0;
        M[N-1] = 0; 
        for (int i = 1; i < N-1; ++i) {
            double h_i = knots[i] - knots[i-1];
            double h_next = knots[i + 1] - knots[i];
            mu[i] = h_i / (h_i + h_next);
            lambda[i] = h_next / (h_i + h_next);
            L[i] = (coefficients[i+1] - coefficients[i]) / (h_next*((h_i + h_next))) - (coefficients[i] - coefficients[i-1]) / (h_i*((h_i + h_next)));
        }

        // 设置三对角矩阵 A 和向量 b
        Eigen::SparseMatrix<double> A(N - 2, N - 2);
        Eigen::VectorXd b(N - 2);

        // 填充三对角矩阵 A 和向量 b
        for (int i = 1; i < N - 3; ++i) {
            A.insert(i, i) = 2;
            A.insert(i, i - 1) = mu[i+1];
            A.insert(i, i + 1) = lambda[i + 1];
            b[i] = 6*L[i+1];
        }
        // 第一个方程的系数
        A.insert(0, 0) = 2; 
        A.insert(0, 1) = lambda[1]; 
        b[0] = 6*L[1] - mu[1]*M[0];

        // 最后一个方程的系数
        A.insert(N-3, N-3) = 2; 
        A.insert(N-3, N-4) = mu[N-2]; 
        b[N-3] = 6*L[N-2] - lambda[N-2]*M[N-1];

        // 使用 Eigen 解线性方程组 A * m = b
        Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
        solver.compute(A);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to decompose the matrix.");
        }
        M0 = solver.solve(b);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to solve the system.");
        }

        for (int i = 1; i < N - 1; ++i){
            M[i] = M0[i-1];
        }

        for (int i = 0; i < N - 1; ++i) {
            if (x >= knots[i] && x < knots[i + 1]) {
                double h_i = knots[i+1] - knots[i]; 
                double K_i = (coefficients[i+1] - coefficients[i])/h_i;
                double m_1 = K_i - (M[i+1]+2*M[i])*h_i/6;
                double m_2 = K_i + (M[i]+2*M[i+1])*h_i/6;
                double sx = coefficients[i] + (x-knots[i])*m_1 + (x-knots[i])*(x-knots[i])*(K_i-m_1)/(knots[i+1]-knots[i]) + (x-knots[i])*(x-knots[i])*(x-knots[i+1])*(m_1+m_2-2*K_i)/((knots[i+1]-knots[i])*(knots[i+1]-knots[i]));
                return sx;
            }
        }
        throw std::out_of_range("x is out of the spline's domain.");
    }

};

// pp形式的完全三次样条(ok)
class ClampedCubicSpline_p : public PPFormSpline {
protected:
    double m1; // 端点处的一阶导数
    double mN; // 端点处的一阶导数
public:
    ClampedCubicSpline_p(const std::vector<double>& knots, const std::vector<double>& coefficients, double m1, double mN)
        : PPFormSpline(knots, coefficients), m1(m1), mN(mN) {}

    // 重写评估函数
    virtual double evaluate(double x) const override{
        assert(knots.size() > 0 && coefficients.size() == knots.size());
    
        int N = knots.size();
        Eigen::VectorXd m(N);
        Eigen::VectorXd m0(N-2);
        Eigen::VectorXd lambda(N);
        Eigen::VectorXd mu(N);
        Eigen::VectorXd K(N);

        // 计算 mu 和 lambda
        m[0] = m1;
        m[N-1] = mN; 
        K[0] = (coefficients[1] - coefficients[0]) / (knots[1] - knots[0]);
        for (int i = 1; i < N-1; ++i) {
            double h_i = knots[i] - knots[i-1];
            double h_next = knots[i + 1] - knots[i];
            mu[i] = h_i / (h_i + h_next);
            lambda[i] = h_next / (h_i + h_next);
            K[i] = (coefficients[i+1] - coefficients[i]) / h_next;
        }

        // 设置三对角矩阵 A 和向量 b
        Eigen::SparseMatrix<double> A(N - 2, N - 2);
        Eigen::VectorXd b(N - 2);

        // 填充三对角矩阵 A 和向量 b
        for (int i = 1; i < N - 3; ++i) {
            A.insert(i, i) = 2;
            A.insert(i, i - 1) = lambda[i+1];
            A.insert(i, i + 1) = mu[i + 1];
            b[i] = 3*mu[i+1]*K[i+1] + 3*lambda[i+1]*K[i];
        }
        // 第一个方程的系数
        A.insert(0, 0) = 2; 
        A.insert(0, 1) = mu[1]; 
        b[0] = 3*mu[1]*K[1] + 3*lambda[1]*K[0] - lambda[1]*m[0];

        // 最后一个方程的系数
        A.insert(N-3, N-3) = 2; 
        A.insert(N-3, N-4) = lambda[N-2]; 
        b[N-3] = 3*mu[N-2]*K[N-2] + 3*lambda[N-2]*K[N-3] - mu[N-2]*m[N-1];

        // 使用 Eigen 解线性方程组 A * m = b
        Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
        solver.compute(A);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to decompose the matrix.");
        }
        m0 = solver.solve(b);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to solve the system.");
        }

        for (int i = 1; i < N - 1; ++i){
            m[i] = m0[i-1];
        }

        for (int i = 0; i < N - 1; ++i) {
            if (x >= knots[i] && x < knots[i + 1]) {
                double sx = coefficients[i] + (x-knots[i])*m[i] + (x-knots[i])*(x-knots[i])*(K[i]-m[i])/(knots[i+1]-knots[i]) + (x-knots[i])*(x-knots[i])*(x-knots[i+1])*(m[i]+m[i+1]-2*K[i])/((knots[i+1]-knots[i])*(knots[i+1]-knots[i]));
                return sx;
            }
        }
        throw std::out_of_range("x is out of the spline's domain.");
    }
};


// 线性B样条函数类(ok)
class LinearSpline_B : public BSpline {
public:
    LinearSpline_B(int degree, const std::vector<double>& knots, const std::vector<double>& coefficients)
        : BSpline(degree, knots, coefficients) {}

    // 重写评估函数
    virtual double evaluate(double x) const override{
        assert(knots.size() > 0 && coefficients.size() == knots.size());
        int N = knots.size();
        for (int i = 0; i < N - 1; ++i) {
            if (x >= knots[i] && x < knots[i + 1]) {
                return coefficients[i+1]*(x-knots[i])/(knots[i+1]-knots[i]) + coefficients[i]*(knots[i+1]-x)/(knots[i+1]-knots[i]);
            }
        }
        throw std::out_of_range("x is out of the spline's domain.");
    }
};

// 二次B样条 (只考虑均匀节点? ok)
class QuadraticSpline_B:  public BSpline {
protected:
    double f1; // 端点处的函数值
    double fN; // 端点处的函数数
public:
    QuadraticSpline_B(int degree, const std::vector<double>& knots, const std::vector<double>& coefficients, double f1, double fN)
        : BSpline(degree, knots, coefficients), f1(f1), fN(fN) {}

    // 计算B样条的二次基函数在给定分段给定节点处的函数值（只考虑均匀节点?）
    double Nip_2(int i, double x0) const{
        double bx;
        double x = x0 + 1.0/2;
        // 这里的B基函数i是在konts[i]左边整点处的
        int N = knots.size(); 
        if(i>0 && i<N-2){
            if(x >= knots[i-1] && x <= knots[i]){
                bx = (x - knots[i-1])*(x - knots[i-1])/((knots[i+1]-knots[i-1])*(knots[i]-knots[i-1]));
            }
            else if(x > knots[i] && x <= knots[i+1]){
                bx = (x - knots[i-1])*(knots[i+1]-x)/((knots[i+1]-knots[i-1])*(knots[i+1]-knots[i]))+(knots[i+2]-x)*(x - knots[i])/((knots[i+2]-knots[i])*(knots[i+1]-knots[i]));
            }
            else if(x > knots[i+1] && x <= knots[i+2]){
                bx = (knots[i+2]-x)*(knots[i+2]-x)/((knots[i+2]-knots[i])*(knots[i+2]-knots[i+1]));
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        else if(i<1){
            if(i==0){
                if(x >= knots[i]-1 && x <= knots[i]){
                    bx = (x - knots[i]+1)*(x - knots[i]+1)/((knots[i+1]-knots[i]+1)*(1));
                }
                else if(x > knots[i] && x <= knots[i+1]){
                    bx = (x - knots[i]+1)*(knots[i+1]-x)/((knots[i+1]-knots[i]+1)*(knots[i+1]-knots[i]))+(knots[i+2]-x)*(x - knots[i])/((knots[i+2]-knots[i])*(knots[i+1]-knots[i]));
                }
                else if(x > knots[i+1] && x <= knots[i+2]){
                    bx = (knots[i+2]-x)*(knots[i+2]-x)/((knots[i+2]-knots[i])*(knots[i+2]-knots[i+1]));
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if(i==-1){
                if(x >= knots[i+1]-2 && x <= knots[i+1]-1){
                    bx = (x - knots[i+1]+2)*(x - knots[i+1]+2)/(2*1);
                }
                else if(x > knots[i+1]-1 && x <= knots[i+1]){
                    bx = (x - knots[i+1]+2)*(knots[i+1]-x)/(2*1)+(knots[i+2]-x)*(x - knots[i+1]+1)/((knots[i+2]-knots[i+1]+1)*(1));
                }
                else if(x > knots[i+1] && x <= knots[i+2]){
                    bx = (knots[i+2]-x)*(knots[i+2]-x)/((knots[i+2]-knots[i+1]+1)*(knots[i+2]-knots[i+1]));
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        else if(i>N-3){
            if(i == N-2){
                if(x >= knots[i-1] && x <= knots[i]){
                    bx = (x - knots[i-1])*(x - knots[i-1])/((knots[i+1]-knots[i-1])*(knots[i]-knots[i-1]));
                }
                else if(x > knots[i] && x <= knots[i+1]){
                    bx = (x - knots[i-1])*(knots[i+1]-x)/((knots[i+1]-knots[i-1])*(knots[i+1]-knots[i]))+(knots[i+1]+1-x)*(x - knots[i])/((knots[i+1]+1-knots[i])*(knots[i+1]-knots[i]));
                }
                else if(x > knots[i+1] && x <= knots[i+1]+1){
                    bx = (knots[i+1]+1-x)*(knots[i+1]+1-x)/((knots[i+1]+1-knots[i])*(1));
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if(i == N-1){
                if(x >= knots[i-1] && x <= knots[i]){
                    bx = (x - knots[i-1])*(x - knots[i-1])/((knots[i]+1-knots[i-1])*(knots[i]-knots[i-1]));
                }
                else if(x > knots[i] && x <= knots[i]+1){
                    bx = (x - knots[i-1])*(knots[i]+1-x)/((knots[i]+1-knots[i-1])*(1))+(knots[i]+2-x)*(x - knots[i])/((2)*(1));
                }
                else if(x > knots[i]+1 && x <= knots[i]+2){
                    bx = (knots[i]+2-x)*(knots[i]+2-x)/((2)*(1));
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if(i == N){
                if(x >= knots[i-1] && x <= knots[i-1]+1){
                    bx = (x - knots[i-1])*(x - knots[i-1])/((2)*(1));
                }
                else if(x > knots[i-1]+1 && x <= knots[i-1]+2){
                    bx = (x - knots[i-1])*(knots[i-1]+2-x)/((2)*(1))+(knots[i-1]+3-x)*(x - knots[i-1]-1)/((2)*(1));
                }
                else if(x > knots[i-1]+2 && x <= knots[i-1]+3){
                    bx = (knots[i-1]+3-x)*(knots[i-1]+3-x)/((2)*(1));
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        else{
            throw std::out_of_range("x is out of the B function's domain.");
        }
        return bx;
    }

    // 重写评估函数
    virtual double evaluate(double x) const override{
        assert(knots.size() > 0 && coefficients.size() == knots.size());
    
        int M = knots.size();
        int N = M+1;
        Eigen::VectorXd a(N+1);
        Eigen::VectorXd a0(N-1);

        // 设置三对角矩阵 A 和向量 b
        Eigen::SparseMatrix<double> A(N-1, N-1);
        Eigen::VectorXd b(N-1);

        // 填充三对角矩阵 A 和向量 b
        for (int i = 1; i < N - 2; ++i) {
            A.insert(i, i) = 6;
            A.insert(i, i - 1) = 1;
            A.insert(i, i + 1) = 1;
            b[i] = 8*coefficients[i];
        }
        // 第一个方程的系数
        A.insert(0, 0) = 5; 
        A.insert(0, 1) = 1; 
        b[0] = 8*coefficients[0] - 2*f1;

        // 最后一个方程的系数
        A.insert(N-2, N-2) = 5; 
        A.insert(N-2, N-3) = 1; 
        b[N-2] = 8*coefficients[N-2] - 2*fN;

        // 使用 Eigen 解线性方程组 A * m = b
        Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;

        // // 观察A的元素值
        // std::cout << "Matrix A (non-zero elements only):" << std::endl;
        // for (int k = 0; k < A.outerSize(); ++k) {
        //     for (Eigen::SparseMatrix<double>::InnerIterator it(A, k); it; ++it) {
        //         std::cout << "A(" << it.row() << "," << it.col() << ") = " << it.value() << std::endl;
        //     }
        // }

        solver.compute(A);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to decompose the matrix.");
        }
        a0 = solver.solve(b);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to solve the system.");
        }

        for (int i = 1; i < N; ++i){
            a[i] = a0[i-1];
        }
        a[0] = 2*f1-a[1];
        a[N] = 2*fN-a[N-1];

        for (int i = 0; i < N - 1; ++i) {
            if (x >= knots[i]- 1.0/2 && x < knots[i]+1.0/2) {
                double sx;
                sx = a[i]*Nip_2(i-1,x) + a[i+1]*Nip_2(i,x) + a[i+2]*Nip_2(i+1,x);
                return sx;
            }
        }
        throw std::out_of_range("x is out of the spline's domain.");
    }
};

// 周期三次B样条(ok)
class PeriodicCubicSpline_B : public BSpline {
public:
    PeriodicCubicSpline_B(int degree, const std::vector<double>& knots, const std::vector<double>& coefficients)
        : BSpline(degree, knots, coefficients) {}
    
    // 计算B样条的三次基函数在给定分段给定x处的函数值
    double Nip_3(int i, double x) const {
        double bx;
        int N = knots.size();
        double t_0_1,t_10,t_21,t_32,t_1_1,t_20,t_2_1,t_31,t_30;
        if(i > 0 && i < N - 3){
            t_0_1 = knots[i] - knots[i-1]; // i>= 1 
            t_10 = knots[i+1] - knots[i];
            t_21 = knots[i+2] - knots[i+1];
            t_32 = knots[i+3] - knots[i+2]; // i<= N-4
            t_1_1 = t_10 + t_0_1;
            t_20 = t_21 + t_10;
            t_2_1 = t_20 + t_0_1;
            t_31 = t_32 + t_21;
            t_30 = t_31 + t_10;
            if(x > knots[i-1] && x <= knots[i]){
                bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
            }
            else if(x > knots[i] && x <= knots[i+1]){
                bx = (x-knots[i-1])/t_2_1 * ((x-knots[i-1])*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+2]-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+3] - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
            }
            else if(x > knots[i+1] && x <= knots[i+2]){
                bx = (x-knots[i-1])*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+3]-x)/t_30*((x-knots[i])*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+3]-x)/(t_31*t_21));
            }
            else if(x > knots[i+2] && x <= knots[i+3]){
                bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        else if (i < 1){
            if(i == -2){
                t_0_1 = 1; // i>= 1 
                t_10 = 1;
                t_21 = 1;
                t_32 = knots[i+3] - knots[i+2]; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i+2] && x < knots[i+3]){
                    bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if(i == -1){
                t_0_1 = 1; // i>= 1 
                t_10 = 1;
                t_21 = knots[i+2] - knots[i+1];
                t_32 = knots[i+3] - knots[i+2]; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i+2] && x < knots[i+3]){
                    bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
                }
                else if(x >= knots[i+1] && x < knots[i+2]){
                    bx = (x-knots[i+1]+2)*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+3]-x)/t_30*((x-knots[i+1]+1)*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+3]-x)/(t_31*t_21));
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if(i == 0){
                t_0_1 = 1; // i>= 1 
                t_10 = knots[i+1] - knots[i];
                t_21 = knots[i+2] - knots[i+1];
                t_32 = knots[i+3] - knots[i+2]; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i+2] && x < knots[i+3]){
                    bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
                }
                else if(x >= knots[i+1] && x < knots[i+2]){
                    bx = (x-knots[i]+1)*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+3]-x)/t_30*((x-knots[i])*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+3]-x)/(t_31*t_21));
                }
                else if(x >= knots[i] && x < knots[i+1]){
                    bx = (x-knots[i]+1)/t_2_1 * ((x-knots[i]+1)*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+2]-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+3] - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        else if (i > N-4){
            if(i == N-1){
                t_0_1 = knots[i] - knots[i-1]; // i>= 1 
                t_10 = 1;
                t_21 = 1;
                t_32 = 1; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x > knots[i-1] && x <= knots[i]){
                    bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if (i == N-2){
                t_0_1 = knots[i] - knots[i-1]; // i>= 1 
                t_10 = knots[i+1] - knots[i];
                t_21 = 1;
                t_32 = 1; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x > knots[i-1] && x <= knots[i]){
                    bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
                }
                else if(x > knots[i] && x <= knots[i+1]){
                    bx = (x-knots[i-1])/t_2_1 * ((x-knots[i-1])*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+1]+1-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+1]+2 - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if (i == N-3){
                t_0_1 = knots[i] - knots[i-1]; // i>= 1 
                t_10 = knots[i+1] - knots[i];
                t_21 = knots[i+2] - knots[i+1];
                t_32 = 1; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x > knots[i-1] && x <= knots[i]){
                    bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
                }
                else if(x > knots[i] && x <= knots[i+1]){
                    bx = (x-knots[i-1])/t_2_1 * ((x-knots[i-1])*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+2]-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+2]+1 - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
                }
                else if(x > knots[i+1] && x <= knots[i+2]){
                    bx = (x-knots[i-1])*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+2]+1-x)/t_30*((x-knots[i])*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+2]+1-x)/(t_31*t_21));
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        return bx; 
    }

    // 重写评估函数
    virtual double evaluate(double x) const override{
        assert(knots.size() > 0 && coefficients.size() == knots.size());
        
        int N = knots.size();
        if ( abs(x - knots[N - 1])< 1e-9 ){
            return coefficients[0];
        }

        Eigen::VectorXd a(N+2);

        // 设置三对角矩阵 A 和向量 b
        Eigen::SparseMatrix<double> A(N+2, N+2);
        Eigen::VectorXd b(N+2);

        // 填充三对角矩阵 A 和向量 b
        for (int i = 1; i < N+1 ; ++i) {
            A.insert(i, i) = Nip_3(i-2,knots[i-1]);
            A.insert(i, i - 1) = Nip_3(i-3,knots[i-1]);
            A.insert(i, i + 1) = Nip_3(i-1,knots[i-1]);
            if(i == N) {
                b[i] = coefficients[0];
            }
            else{
                b[i] = coefficients[i-1];
            }
        }
        // 第一个方程的系数
        double c_1,c0,c1;
        double d_1,d0,d1;
        c_1 = -(knots[1]-knots[0])/((knots[1]-knots[0]+2)*(knots[1]-knots[0]+1));
        c1 = 1/((knots[2]-knots[0]+1)*(knots[1]-knots[0]+1));
        c0 = -c_1 -c1;
        d_1 = -1/((knots[N-1]-knots[N-3]+1)*(knots[N-1]-knots[N-2]+1));
        d1 = (knots[N-1]-knots[N-2])/((knots[N-1]-knots[N-2]+2)*(knots[N-1]-knots[N-2]+1));
        d0 = -d_1-d1;
        A.insert(0, 0) = c_1; 
        A.insert(0, 1) = c0; 
        A.insert(0, 2) = c1;
        A.insert(0, N-1) = -d_1; 
        A.insert(0, N) = -d0; 
        A.insert(0, N+1) = -d1;  

        // 最后一个方程的系数
        double e_1,e0,e1;
        double f_1,f0,f1;
        e_1 = 1/((knots[1]-knots[0]+2)*(knots[1]-knots[0]+1));
        e1 = 1/((knots[2]-knots[0]+1)*(knots[1]-knots[0]+1));
        e0 = -e_1 -e1;
        f_1 = 1/((knots[N-1]-knots[N-3]+1)*(knots[N-1]-knots[N-2]+1));
        f1 = 1/((knots[N-1]-knots[N-2]+2)*(knots[N-1]-knots[N-2]+1));
        f0 = -f_1 -f1;
        A.insert(N+1, 0) = e_1; 
        A.insert(N+1, 1) = e0; 
        A.insert(N+1, 2) = e1;
        A.insert(N+1, N-1) = -f_1; 
        A.insert(N+1, N) = -f0; 
        A.insert(N+1, N+1) = -f1; 

        // 使用 Eigen 解线性方程组 A * m = b
        Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;

        // // 观察A的元素值
        // std::cout << "Matrix A (non-zero elements only):" << std::endl;
        // for (int k = 0; k < A.outerSize(); ++k) {
        //     for (Eigen::SparseMatrix<double>::InnerIterator it(A, k); it; ++it) {
        //         std::cout << "A(" << it.row() << "," << it.col() << ") = " << it.value() << std::endl;
        //     }
        // }

        solver.compute(A);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to decompose the matrix.");
        }
        a = solver.solve(b);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to solve the system.");
        }

        for (int i = 0; i < N - 1; ++i) {
            if (x >= knots[i] && x < knots[i + 1]) {
                double sx;
                if(x == knots[i]) { sx = coefficients[i];}
                else {
                    sx = a[i]*Nip_3(i-2,x) + a[i+1]*Nip_3(i-1,x) + a[i+2]*Nip_3(i,x) + a[i+3]*Nip_3(i+1,x);
                }
                return sx;
            }
        }
        throw std::out_of_range("x is out of the spline's domain.");
    }

};

// 自然三次B样条(ok)
class NaturalSpline_B : public BSpline {
public:
    NaturalSpline_B(int degree, const std::vector<double>& knots, const std::vector<double>& coefficients)
        : BSpline(degree, knots, coefficients) {}
    // 计算B样条的三次基函数在给定分段给定x处的函数值
    double Nip_3(int i, double x) const {
        double bx;
        int N = knots.size();
        double t_0_1,t_10,t_21,t_32,t_1_1,t_20,t_2_1,t_31,t_30;
        if(i > 0 && i < N - 3){
            t_0_1 = knots[i] - knots[i-1]; // i>= 1 
            t_10 = knots[i+1] - knots[i];
            t_21 = knots[i+2] - knots[i+1];
            t_32 = knots[i+3] - knots[i+2]; // i<= N-4
            t_1_1 = t_10 + t_0_1;
            t_20 = t_21 + t_10;
            t_2_1 = t_20 + t_0_1;
            t_31 = t_32 + t_21;
            t_30 = t_31 + t_10;
            if(x > knots[i-1] && x <= knots[i]){
                bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
            }
            else if(x > knots[i] && x <= knots[i+1]){
                bx = (x-knots[i-1])/t_2_1 * ((x-knots[i-1])*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+2]-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+3] - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
            }
            else if(x > knots[i+1] && x <= knots[i+2]){
                bx = (x-knots[i-1])*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+3]-x)/t_30*((x-knots[i])*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+3]-x)/(t_31*t_21));
            }
            else if(x > knots[i+2] && x <= knots[i+3]){
                bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        else if (i < 1){
            if(i == -2){
                t_0_1 = 1; // i>= 1 
                t_10 = 1;
                t_21 = 1;
                t_32 = knots[i+3] - knots[i+2]; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i+2] && x < knots[i+3]){
                    bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if(i == -1){
                t_0_1 = 1; // i>= 1 
                t_10 = 1;
                t_21 = knots[i+2] - knots[i+1];
                t_32 = knots[i+3] - knots[i+2]; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i+2] && x < knots[i+3]){
                    bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
                }
                else if(x >= knots[i+1] && x < knots[i+2]){
                    bx = (x-knots[i+1]+2)*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+3]-x)/t_30*((x-knots[i+1]+1)*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+3]-x)/(t_31*t_21));
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if(i == 0){
                t_0_1 = 1; // i>= 1 
                t_10 = knots[i+1] - knots[i];
                t_21 = knots[i+2] - knots[i+1];
                t_32 = knots[i+3] - knots[i+2]; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i+2] && x < knots[i+3]){
                    bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
                }
                else if(x >= knots[i+1] && x < knots[i+2]){
                    bx = (x-knots[i]+1)*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+3]-x)/t_30*((x-knots[i])*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+3]-x)/(t_31*t_21));
                }
                else if(x >= knots[i] && x < knots[i+1]){
                    bx = (x-knots[i]+1)/t_2_1 * ((x-knots[i]+1)*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+2]-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+3] - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        else if (i > N-4){
            if(i == N-1){
                t_0_1 = knots[i] - knots[i-1]; // i>= 1 
                t_10 = 1;
                t_21 = 1;
                t_32 = 1; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x > knots[i-1] && x <= knots[i]){
                    bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if (i == N-2){
                t_0_1 = knots[i] - knots[i-1]; // i>= 1 
                t_10 = knots[i+1] - knots[i];
                t_21 = 1;
                t_32 = 1; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x > knots[i-1] && x <= knots[i]){
                    bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
                }
                else if(x > knots[i] && x <= knots[i+1]){
                    bx = (x-knots[i-1])/t_2_1 * ((x-knots[i-1])*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+1]+1-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+1]+2 - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if (i == N-3){
                t_0_1 = knots[i] - knots[i-1]; // i>= 1 
                t_10 = knots[i+1] - knots[i];
                t_21 = knots[i+2] - knots[i+1];
                t_32 = 1; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x > knots[i-1] && x <= knots[i]){
                    bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
                }
                else if(x > knots[i] && x <= knots[i+1]){
                    bx = (x-knots[i-1])/t_2_1 * ((x-knots[i-1])*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+2]-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+2]+1 - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
                }
                else if(x > knots[i+1] && x <= knots[i+2]){
                    bx = (x-knots[i-1])*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+2]+1-x)/t_30*((x-knots[i])*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+2]+1-x)/(t_31*t_21));
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        return bx; 
    }

    // 重写评估函数
    virtual double evaluate(double x) const override{
        assert(knots.size() > 0 && coefficients.size() == knots.size());
    
        int N = knots.size();
        Eigen::VectorXd a(N+2);
        Eigen::VectorXd a0(N);

        // 设置三对角矩阵 A 和向量 b
        Eigen::SparseMatrix<double> A(N, N);
        Eigen::VectorXd b(N);

        // 填充三对角矩阵 A 和向量 b
        for (int i = 1; i < N - 1; ++i) {
            A.insert(i, i) = Nip_3(i-1,knots[i]);
            A.insert(i, i - 1) = Nip_3(i-2,knots[i]);
            A.insert(i, i + 1) = Nip_3(i,knots[i]);
            b[i] = coefficients[i];
        }
        // 第一个方程的系数
        double c0,c1;
        c1 = -(knots[1]-knots[0]+2)/(knots[2]-knots[0]+1);
        c0 = 1-c1;
        A.insert(0, 0) = Nip_3(-1,knots[0]) + Nip_3(-2,knots[0])*c0; 
        A.insert(0, 1) = Nip_3(0,knots[0]) + Nip_3(-2,knots[0])*c1; 
        b[0] = coefficients[0];

        // 最后一个方程的系数
        double d0,d1;
        d1 = -(knots[N-1]-knots[N-2]+2)/(knots[N-1]-knots[N-3]+1);
        d0 = 1-d1;
        A.insert(N-1, N-1) = Nip_3(N-2,knots[N-1])+ Nip_3(N-1,knots[N-1])*d0; 
        A.insert(N-1, N-2) = Nip_3(N-3,knots[N-1])+ Nip_3(N-1,knots[N-1])*d1; 
        b[N-1] = coefficients[N-1];

        // 使用 Eigen 解线性方程组 A * m = b
        Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;

        // // 观察A的元素值
        // std::cout << "Matrix A (non-zero elements only):" << std::endl;
        // for (int k = 0; k < A.outerSize(); ++k) {
        //     for (Eigen::SparseMatrix<double>::InnerIterator it(A, k); it; ++it) {
        //         std::cout << "A(" << it.row() << "," << it.col() << ") = " << it.value() << std::endl;
        //     }
        // }

        solver.compute(A);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to decompose the matrix.");
        }
        a0 = solver.solve(b);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to solve the system.");
        }

        for (int i = 1; i < N + 1; ++i){
            a[i] = a0[i-1];
        }
        a[0] = c0*a[1] + c1*a[2];
        a[N+1] = d0*a[N] + d1*a[N-1];

        for (int i = 0; i < N - 1; ++i) {
            if (x >= knots[i] && x < knots[i + 1]) {
                double sx;
                if(x == knots[i]) { sx = coefficients[i];}
                else {
                    sx = a[i]*Nip_3(i-2,x) + a[i+1]*Nip_3(i-1,x) + a[i+2]*Nip_3(i,x) + a[i+3]*Nip_3(i+1,x);
                }
                return sx;
            }
        }
        throw std::out_of_range("x is out of the spline's domain.");
    }

};

// 完全三次B样条(ok)
class ClampedCubicSpline_B : public BSpline {
protected:
    double m1; // 端点处的一阶导数
    double mN; // 端点处的一阶导数
public:
    ClampedCubicSpline_B(int degree, const std::vector<double>& knots, const std::vector<double>& coefficients, double m1, double mN)
        : BSpline(degree, knots, coefficients), m1(m1), mN(mN) {}

    // 计算B样条的三次基函数在给定分段给定x处的函数值
    double Nip_3(int i, double x) const {
        double bx;
        int N = knots.size();
        double t_0_1,t_10,t_21,t_32,t_1_1,t_20,t_2_1,t_31,t_30;
        if(i > 0 && i < N - 3){
            t_0_1 = knots[i] - knots[i-1]; // i>= 1 
            t_10 = knots[i+1] - knots[i];
            t_21 = knots[i+2] - knots[i+1];
            t_32 = knots[i+3] - knots[i+2]; // i<= N-4
            t_1_1 = t_10 + t_0_1;
            t_20 = t_21 + t_10;
            t_2_1 = t_20 + t_0_1;
            t_31 = t_32 + t_21;
            t_30 = t_31 + t_10;
            if(x >= knots[i-1] && x <= knots[i]){
                bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
            }
            else if(x > knots[i] && x <= knots[i+1]){
                bx = (x-knots[i-1])/t_2_1 * ((x-knots[i-1])*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+2]-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+3] - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
            }
            else if(x > knots[i+1] && x <= knots[i+2]){
                bx = (x-knots[i-1])*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+3]-x)/t_30*((x-knots[i])*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+3]-x)/(t_31*t_21));
            }
            else if(x > knots[i+2] && x <= knots[i+3]){
                bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        else if (i < 1){
            if(i == -2){
                t_0_1 = 1; // i>= 1 
                t_10 = 1;
                t_21 = 1;
                t_32 = knots[i+3] - knots[i+2]; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i+2] && x < knots[i+3]){
                    bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if(i == -1){
                t_0_1 = 1; // i>= 1 
                t_10 = 1;
                t_21 = knots[i+2] - knots[i+1];
                t_32 = knots[i+3] - knots[i+2]; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i+2] && x < knots[i+3]){
                    bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
                }
                else if(x >= knots[i+1] && x < knots[i+2]){
                    bx = (x-knots[i+1]+2)*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+3]-x)/t_30*((x-knots[i+1]+1)*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+3]-x)/(t_31*t_21));
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if(i == 0){
                t_0_1 = 1; // i>= 1 
                t_10 = knots[i+1] - knots[i];
                t_21 = knots[i+2] - knots[i+1];
                t_32 = knots[i+3] - knots[i+2]; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i+2] && x < knots[i+3]){
                    bx = (knots[i+3]-x)*(knots[i+3]-x)*(knots[i+3]-x)/(t_30*t_31*t_32);
                }
                else if(x >= knots[i+1] && x < knots[i+2]){
                    bx = (x-knots[i]+1)*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+3]-x)/t_30*((x-knots[i])*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+3]-x)/(t_31*t_21));
                }
                else if(x >= knots[i] && x < knots[i+1]){
                    bx = (x-knots[i]+1)/t_2_1 * ((x-knots[i]+1)*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+2]-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+3] - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        else if (i > N-4){
            if(i == N-1){
                t_0_1 = knots[i] - knots[i-1]; // i>= 1 
                t_10 = 1;
                t_21 = 1;
                t_32 = 1; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i-1] && x <= knots[i]){
                    bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if (i == N-2){
                t_0_1 = knots[i] - knots[i-1]; // i>= 1 
                t_10 = knots[i+1] - knots[i];
                t_21 = 1;
                t_32 = 1; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i-1] && x <= knots[i]){
                    bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
                }
                else if(x > knots[i] && x <= knots[i+1]){
                    bx = (x-knots[i-1])/t_2_1 * ((x-knots[i-1])*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+1]+1-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+1]+2 - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else if (i == N-3){
                t_0_1 = knots[i] - knots[i-1]; // i>= 1 
                t_10 = knots[i+1] - knots[i];
                t_21 = knots[i+2] - knots[i+1];
                t_32 = 1; // i<= N-4
                t_1_1 = t_10 + t_0_1;
                t_20 = t_21 + t_10;
                t_2_1 = t_20 + t_0_1;
                t_31 = t_32 + t_21;
                t_30 = t_31 + t_10;
                if(x >= knots[i-1] && x <= knots[i]){
                    bx = (x - knots[i-1])*(x - knots[i-1])*(x - knots[i-1])/(t_2_1*t_1_1*t_0_1);
                }
                else if(x > knots[i] && x <= knots[i+1]){
                    bx = (x-knots[i-1])/t_2_1 * ((x-knots[i-1])*(knots[i+1]-x)/(t_1_1*t_10)+(knots[i+2]-x)*(x-knots[i])/(t_20*t_10)) + (knots[i+2]+1 - x)*(x-knots[i])*(x-knots[i])/(t_30*t_20*t_10);
                }
                else if(x > knots[i+1] && x <= knots[i+2]){
                    bx = (x-knots[i-1])*(knots[i+2]-x)*(knots[i+2]-x)/(t_2_1*t_20*t_21) + (knots[i+2]+1-x)/t_30*((x-knots[i])*(knots[i+2]-x)/(t_20*t_21)+(x-knots[i+1])*(knots[i+2]+1-x)/(t_31*t_21));
                }
                else{
                    throw std::out_of_range("x is out of the B function's domain.");
                }
            }
            else{
                throw std::out_of_range("x is out of the B function's domain.");
            }
        }
        return bx; 
    }
    
    // 重写评估函数
    virtual double evaluate(double x) const override{
        assert(knots.size() > 0 && coefficients.size() == knots.size());
    
        int N = knots.size();
        Eigen::VectorXd a(N+2);
        Eigen::VectorXd a0(N);

        // 设置三对角矩阵 A 和向量 b
        Eigen::SparseMatrix<double> A(N, N);
        Eigen::VectorXd b(N);

        // 填充三对角矩阵 A 和向量 b
        for (int i = 1; i < N - 1; ++i) {
            A.insert(i, i) = Nip_3(i-1,knots[i]);
            A.insert(i, i - 1) = Nip_3(i-2,knots[i]);
            A.insert(i, i + 1) = Nip_3(i,knots[i]);
            b[i] = coefficients[i];
        }
        // 第一个方程的系数
        double c0,c1,c2;
        c1 = (knots[1]-knots[0]+2)/((knots[1]-knots[0])*(knots[2]-knots[0]+1));
        c0 = 1-c1;
        c2 = ((knots[1]-knots[0]+2)*(knots[1]-knots[0]+1))/(3*(knots[1]-knots[0]));
        A.insert(0, 0) = Nip_3(-1,knots[0]) + Nip_3(-2,knots[0])*c0; 
        A.insert(0, 1) = Nip_3(0,knots[0]) + Nip_3(-2,knots[0])*c1; 
        b[0] = coefficients[0] + Nip_3(-2,knots[0])*c2*m1;

        // 最后一个方程的系数
        double d0,d1,d2;
        d1 = (knots[N-1]-knots[N-2]+2)/((knots[N-1]-knots[N-3]+1)*(knots[N-1]-knots[N-2]));
        d0 = 1-d1;
        d2 = ((knots[N-1]-knots[N-2]+2)*(knots[N-1]-knots[N-2]+1))/(3*(knots[N-1]-knots[N-2]));
        A.insert(N-1, N-1) = Nip_3(N-2,knots[N-1])+ Nip_3(N-1,knots[N-1])*d0; 
        A.insert(N-1, N-2) = Nip_3(N-3,knots[N-1])+ Nip_3(N-1,knots[N-1])*d1; 
        b[N-1] = coefficients[N-1] - Nip_3(N-1,knots[N-1])*d2*mN;

        // 使用 Eigen 解线性方程组 A * m = b
        Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;

        // // 观察A的元素值
        // std::cout << "Matrix A (non-zero elements only):" << std::endl;
        // for (int k = 0; k < A.outerSize(); ++k) {
        //     for (Eigen::SparseMatrix<double>::InnerIterator it(A, k); it; ++it) {
        //         std::cout << "A(" << it.row() << "," << it.col() << ") = " << it.value() << std::endl;
        //     }
        // }

        solver.compute(A);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to decompose the matrix.");
        }
        a0 = solver.solve(b);
        if (solver.info() != Eigen::Success) {
            throw std::runtime_error("Failed to solve the system.");
        }

        for (int i = 1; i < N + 1; ++i){
            a[i] = a0[i-1];
        }
        a[0] = c0*a[1] + c1*a[2] - c2*m1 ;
        a[N+1] = d0*a[N] + d1*a[N-1]+ d2*mN;

        for (int i = 0; i < N - 1; ++i) {
            if (x >= knots[i] && x < knots[i + 1]) {
                double sx;
                
                    sx = a[i]*Nip_3(i-2,x) + a[i+1]*Nip_3(i-1,x) + a[i+2]*Nip_3(i,x) + a[i+3]*Nip_3(i+1,x);
                
                return sx;
            }
        }
        throw std::out_of_range("x is out of the spline's domain.");
    }

};

// 任意阶B样条函数(ok)
class ArbitraryOrderBSpline : public BSpline {
public:
    ArbitraryOrderBSpline(int degree, const std::vector<double>& knots, const std::vector<double>& coefficients)
        : BSpline(degree, knots, coefficients) {}

    // 重写评估函数
    virtual double evaluate(double x) const override {
        // 这里的coefficients不再是节点函数值，而是从2-n到N的基函数的系数
        assert(knots.size() > 0 && coefficients.size() == degree + knots.size() - 1);

        int N = knots.size(); // 计算节点数量
        int k = N-1;

        // 找出x所在区间
        for(int j = 0; j< N-1; j++){
            if(x >= knots[j] && x < knots[j + 1]) {
                k = j;
                break;
            }
        }
        if( k == N-1 && x!=knots[N-1]){
            throw std::out_of_range("x is out of the spline function's domain.");
        }

        double sum = 0; // 初始化求和
        for (int i = std::max( 1-degree, k - degree + 1); i < std::min(k+1,N-1) + 1 ; ++i) {
            sum += coefficients[i-(1 - degree)] * Bspline_value(i, degree, x);
        }
        return sum;
    }

protected:
    // 计算B样条基函数值的递归函数
    double Bspline_value(int i, int n, double x) const {
        int N = knots.size();
        // 很蠢的延拓法，把超出边界的节点设置为+-1
        if (n == 0) {
            if(i < 0) {
                return (x > knots[0]+i-1 && x <= knots[0]+i) ? 1.0 : 0.0;
            } else if(i < 1) {
                return (x > knots[0]+i-1 && x <= knots[i]) ? 1.0 : 0.0;
            } else if(i < N){
                return (x > knots[i-1] && x <= knots[i]) ? 1.0 : 0.0; // 初始版本
            } else if(i < N+1){
                return (x > knots[i-1] && x <= knots[N-1]+i-(N-1)) ? 1.0 : 0.0;
            } else {
                return (x > knots[N-1]+i-1-(N-1) && x <= knots[N-1]+i-(N-1)) ? 1.0 : 0.0;
            }
        } else {
            double a,b;
            if(i<0){
                a = (x - knots[0]-(i-1)) / (knots[i + n - 1] - knots[0]-(i-1)); 
                b = (knots[i + n] - x) / (knots[i + n] - knots[0]-i);
            } else if(i<1){
                a = (x - knots[0]-(i-1)) / (knots[i + n - 1] - knots[0]-(i-1));
                b = (knots[i + n] - x) / (knots[i + n] - knots[i]);
            } else if(i<N-n){
                a = (x - knots[i-1]) / (knots[i + n - 1] - knots[i - 1]); //初始版本
                b = (knots[i + n] - x) / (knots[i + n] - knots[i]);
            } else if(i<N-n+1){
                a = (x - knots[i-1]) / (knots[i + n - 1] - knots[i - 1]);
                b = (knots[N-1]+i+n-(N-1) - x) / (knots[N-1]+i+n-(N-1) - knots[i]);
            } else {
                a = (x - knots[i-1]) / (knots[N-1]+i+n-1-(N-1) - knots[i - 1]);
                b = (knots[N-1]+i+n-(N-1) - x) / (knots[N-1]+i+n-(N-1) - knots[i]);
            }
            return a * Bspline_value(i, n - 1, x) + b * Bspline_value(i + 1, n - 1, x);
        }
    }
};

#endif // SPLINE_FUNCTIONS_H